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Certain  commercial  software  products  are  identified  in  this  paper  in  order  to  adequately 
specify  the  computational  procedures.  Such  identification  does  not  imply 
recommendation  or  endorsement  by  the  National  Institute  of  Standards  and  Technology, 
nor  does  it  imply  that  the  software  products  identified  are  necessarily  the  best  available 
for  the  purpose. 


Abstract 


The  precision  of  manufacturing  using  machine  tools  depends  on  the  accuracy  of  the 
relative  position  of  the  cutting  tool  with  respect  to  the  workpiece.  Kinematic  modeling  of 
machine  tools  is  used  to  describe  this  relative  position.  This  motion  can  be  modeled  by 
homogeneous  coordinate  transformation  matrices  composed  of  both  rotational  elements 
as  well  as  positional  offset  elements  of  the  associated  coordinates.  The  rotation  and 
translation  components  of  the  homogeneous  transformations  are  considered  to  be 
functions  of  nominal  tool  position  and  machine  temperatures.  In  general  these  functions 
are  low  order  polynomials  in  terms  of  position  and  temperatures.  The  coefficients  are 
usually  calculated  with  least  squares  curve  fitting  techniques.  The  data  for  these  fits  are 
obtained  by  measuring  actual  coordinate  positions  and  temperatures  on  the  machine  tool 
based  upon  desired  programmed  nominal  coordinates.  The  process  of  measuring  and 
modeling  the  errors  of  machine  axis  positions  as  functions  of  nominal  positions  and 
temperatures  is  referred  to  as  machine  tool  characterization.  The  geometric-thermal 
models  developed  through  machine  tool  characterization  may  not  fully  predict  the  errors 
encountered  by  a machine  tool  during  machining.  The  data  from  machine 
characterization  usually  provides  the  structure  to  derive  the  basic  form  of  the  equations 
used  to  model  the  various  error  components  used  in  the  homogeneous  matrices.  This 
process  of  model  updating  involves  determining  the  residual  systematic  errors  of  the 
machine  tool  and  applying  an  algorithm  to  update  the  geometric-thermal  model 
coefficients.  The  updating  algorithm  described  in  this  report  begins  with  adding 
perturbation  terms  to  the  characterization  coefficients  of  the  geometric-thermal  model. 
These  coefficients  are  estimated  by  an  “inverse”  process,  using  residual  systematic  errors, 
determined  from  part  measurements  on  a coordinate  measuring  machine.  The  main  tool 
used  in  identifying  the  perturbation  terms  is  called  a generalized  or  pseudo  inverse 
matrix.  This  matrix  is  applied  to  the  residual  error  vector  to  obtain  a “best”  approximate 
solution  to  the  least  squares  problem. 

Key  Words:  generalized  inverse;  geometric-thermal  model;  least  squares;  machine  tool 
model;  post-process  analysis;  pseudo  inverse 
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0.0  Introduction 


The  precision  of  manufacturing  using  machine  tools  depends  on  the  accuracy  of  the 
relative  position  of  the  cutting  tool  with  respect  to  the  workpiece.  Kinematic  modeling  of 
machine  tools  is  used  to  describe  this  relative  position.  In  this  paper  a machine  tool  is 
assumed  to  be  decomposable  into  basic  rigid  components.  The  motions  of  these 
components  can  be  related  to  each  other  by  means  of  relative  motions  of  coordinate 
systems  associated  with  each  component.  The  machine  tool  is  assumed  to  have  a 
reference  axis  system  (arbitrarily  assigned  by  the  user).  From  this  axis  system  two  chains 
of  linked  coordinates  emanate.  One  connects  the  reference  axis  system  with  the  tool 
point  and  the  other  with  the  workpiece.  The  motion  of  a coordinate  system  relative  to 
another  can  be  modeled  by  a homogeneous  coordinate  transformation  matrix.  This 
matrix  is  composed  of  both  rotational  elements  as  well  as  positional  offset  elements  of 
the  associated  coordinates.  The  significance  of  homogeneous  matrices  is  that  a 
coordinate  system  defining  the  motion  of  a machine  component  can  be  linked  to  the 
reference  coordinate  system  by  a sequence  of  matrix  multiplications.  The  errors  in  the 
location  of  the  tool  with  respect  to  the  workpiece  are  related  to  the  relative  errors  in  the 
locations  of  the  coordinate  frames  in  the  axis  chains.  The  rotation  and  translation 
components  of  the  homogeneous  matrices  defining  the  relative  motions  of  the  linked 
coordinate  systems  are  considered  to  be  functions  of  nominal  tool  position  and  machine 
temperatures.  In  general  these  functions  are  low  order  polynomials  in  terms  of  position 
and  temperatures.  The  coefficients  are  usually  calculated  with  least  squares  curve  fitting 
techniques.  The  data  for  these  fits  are  obtained  by  measuring  actual  coordinate  positions 
and  temperatures  on  the  machine  tool  based  upon  desired  (CNC  programmed)  nominal 
coordinates.  Modeling  of  machine  tools  in  this  manner  has  been  used  by  many 
researchers,  see  for  example  references  [1]  to  [1 1]. 

The  process  of  measuring  and  modeling  the  errors  of  machine  axis  positions  as  functions 
of  nominal  positions  and  temperatures  is  referred  to  as  machine  tool  characterization. 

One  of  the  principal  applications  of  kinematic  modeling  and  machine  tool 
characterization  is  to  develop  a geometric-thermal  model  that  relates  the  machine  tool 
error  to  the  geometric  imperfection  and  temperature  variation  in  the  machine  tool.  This 
model  can  also  be  used  for  real-time  correction  of  machine  tool  errors  as  functions  of 
nominal  position  and  machine  tool  temperature  profile.  References  [3]  and  [12]  through 
[15]  describe  several  machine  tool  characterization  experiments. 

The  geometric-thermal  models  developed  through  machine  tool  characterization  may  not 
fully  predict  the  errors  encountered  by  a machine  tool  during  machining.  This  is  because 
the  characterization  is  usually  based  on  measurements  that  can  be  taken  only  when  the 
machine  tool  is  not  operating.  Many  factors,  such  as  coolant  spray,  load-induced 
deformation,  and  long  term  machine  wear,  combine  to  make  such  models  ineffective. 

The  data  from  machine  characterization  usually  provides  the  structure  to  derive  the  basic 
form  of  the  equations  used  to  model  the  various  error  components  used  in  the 
homogeneous  matrices.  The  models  developed  through  machine  characterization  can  be 
considered  as  first  order  models  that  require  coefficient  adjustment  and  periodic 
coefficient  updating  in  order  to  capture  slowly  changing  process  conditions  in  a 
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production  environment.  This  process  of  model  updating  involves  determining  the 
residual  systematic  errors  of  the  machine  tool  and  applying  an  algorithm  to  update  the 
geometric-thermal  model  coefficients. 

The  updating  algorithm  described  in  this  report  begins  with  adding  perturbation  terms  to 
the  characterization  coefficients  of  the  geometric-thermal  model.  These  coefficients  are 
estimated  by  an  “inverse”  process,  using  residual  systematic  errors,  determined  from  part 
measurements  on  a coordinate  measuring  machine  (CMM).  They  are  then  combined 
with  the  existing  model  coefficients  to  compute  the  updated  model  coefficients.  The 
“inverse”  process  used  is  a form  of  least-squares  estimation  that  takes  into  account  the 
possibility  that  the  equations  that  arise  in  the  usual  least  squares  process  may  not  be 
uniquely  solvable  and  standard  least  squares  procedures  may  fail. 

A discussion  of  one  implementation  of  Post-Process  analysis  can  be  found  in  [20] 
through  [23]  although  few  details  of  the  computational  inversion  process  are  given.  That 
implementation  relies  on  inversion  of  data  acquired  on  individual  features.  The  current 
algorithm  significantly  differs  from  that  implementation  in  that  it  uses  data  obtained  from 
measurements  of  an  entire  part  to  determine  residual  errors  rather  than  restricting  the 
inversion  process  to  a feature  at  a time.  Using  data  from  an  entire  part  provides  an 
averaged  set  of  coefficients  applicable  to  all  features  on  the  part.  Inverting,  based  on 
fitting  separate  feature  measurements,  adjusts  the  same  coefficients  in  the  geometric- 
thermal  model  differently  and  does  not  provide  a single  consistent  adjustment  of  all 
parameters  that  apply  to  the  entire  part.  Once  all  parameters  are  uniformly  adjusted, 
errors  on  individual  features  can  easily  be  determined. 

The  main  tool  used  in  identifying  the  perturbation  terms  is  called  a generalized  or  pseudo 
inverse  matrix.  This  matrix  is  applied  to  the  residual  error  vector  to  obtain  a “best” 
approximate  solution  to  the  least  squares  problem.  The  term  “best”  will  be  described  in 
Appendix  E.  The  essential  tool  in  forming  the  generalized  inverse  is  a matrix  theorem 
which  proves  that  any  given  matrix  can  be  decomposed  into  an  associated  diagonal 
matrix,  some  of  whose  elements  might  be  zero.  The  diagonal  elements  are  called 
singular  values  and  indicate  how  many  perturbation  coefficients  are  significant.  A full 
discussion  of  the  numerical  aspects  related  to  the  singular  value  decomposition  and 
pseudo  inverse  generation  are  given  in  Appendix  E and  in  reference  [24], 

This  report  describes  the  algorithmic  details  of  the  post-process  analysis  procedure  and 
presents  an  application  to  one  part  produced  on  a turning  center.  The  authors  are  aware 
that  the  application  of  the  method  to  one  part  and,  in  fact,  a reasonably  simple  part  does 
not  constitute  an  adequate  test  of  the  procedure.  However,  budgetary  constraints  at  the 
time  of  the  study  limited  the  current  work.  The  authors,  however,  felt  that  it  was 
worthwhile  publishing  the  current  results  as  a stimulus  for  further  research.  Future  testing 
of  the  algorithm  is  planned  as  budgets  allow.  These  new  tests  will  highlight  various  parts 
designed  to  utilize  the  extent  of  the  workspace  and  temperature  ranges  experienced  by  the 
turning  center.  The  object  of  the  tests  is  to  determine  whether  there  is  a fixed  set  of 
coefficients  for  the  geometric-thermal  model  that  apply  to  the  entire  workspace  or 
whether  the  workspace  needs  to  be  divided  into  sub-areas  with  separate  sets  of 
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coefficients  applicable  to  each  sub-area.  They  will  also  study  the  robustness  of  the 
compensated  parameters  to  perturbations  of  the  data  points  upon  which  the  model 
updates  were  based. 


The  report  is  divided  into  six  sections.  Section  2 presents  the  error  equations  for  the 
machine  tool  (in  this  report  a turning  center  was  used).  Section  3 describes  the 
development  of  the  perturbation  equations.  Section  3 describes  the  formation  of  the 
system  matrix.  Section  4 describes  the  application  of  the  method  to  estimate  the 
perturbed  machine  tool  model  coefficients  by  using  a machined  part  on  the  turning 
center.  Section  5 concludes  with  a discussion  of  the  results.  Section  6 includes  the 
references. 


7 


1.0  Kinematic  Model  for  a Turning  Center 


The  turning  center  that  is  modeled  in  this  section  is  shown  in  Figure  1.  Although  40 
thermocouples  have  been  identified  numerically  from  0 to  39  (see  Appendix  A.l,  Table 
Al),  only  36  were  actually  used  to  measure  temperatures  since  3 through  6 were  not 
attached  to  the  machine  tool  during  testing.  This  figure  also  shows  the  locations  of 
various  thermocouples  used  to  model  the  individual  component  errors  that  are  included  in 
the  complete  kinematic  model. 


Figure  1:  Turning  Center  with  Thermocouple  Locations 


The  formal  development  of  the  geometric-thermal  error  equations  from  the  kinematic 
model  is  given  in  Appendix  A.l  and  only  the  results  are  given  below.  The  definitions  of 
the  terms  in  the  model  are  also  given  in  Appendix  A.L  The  methods  used  to  develop 
these  equations  are  well  known  and  can  be  found  in  [3-5].  The  coordinate  system  for  the 
turning  center  is  shown  in  Figure  1. 
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The  modeled  displacement  errors  of  the  tool  from  the  workpiece  are  given  by: 


Ex=ey(s.T)*Zs+{ey{z,T)  + ey(x,T)  + ey(t,T)-ey(s,T)\*Z. 

+ [ey(z,T)~ £y(s,T)]*Z2  - ey(s,T)*(Z , + z)  + Sx(z,T)  + Sx(x,T) 
- Sx  (w,  T)  + Sx(t,T)~  Sx  (s,T)  + Sx  (c,  T) 


Ey=-£x(s,T)*Z5+[-£x(z,T)-£x(x,T)-£x(t,T)  + £x(s,T)]*Z4 
+ [-£x(z,T)  + £x(s,T)]*Z2+£x(s,T)*(Zi+z)  + [£z(z,T)  + £z(x,T)  + £:(t,T)]*X4  (1) 

+ £z(z,T)*(X2+x)  + Sy(z,T)  + Sy(x,T)-5y(w,T)  + Sy(t,T)-Sy(s,T)  + Sy(c,T ) 

Ez  =-£y(s,T)*X5+[-£y(z,T)-£y(x,T)-£y(t,T)  + £y(s,T)]*XA 
+ [-£y(z,T)  + £y(s,T)]*(X2+x)  + £y(s,T)*X]  +S.(z,T)  + Sz  (x,T) 

- S.  (w,T)  + Sz  (t,T)-Sz(s,T)  + S.  (c,  T) 

These  equations  give  the  generic  error  formulas  for  a turning  center  with  its  absolute 
reference  axis  system  and  machine  coordinate  system  assumed  aligned  at  “cold  start”  in 
the  “home”  position.  The  constants,  X,,  Z,,  X2,  Z2,  X4,  Z4,  X5,  Z5 , are  relative  offsets 
of  the  machine  component  coordinate  systems  from  the  absolute  reference  axis  system. 
Note  that  X3,  Z3  are  found  to  be  non-significant  and  thus  do  not  appear  in  equation  (1). 

The  significant  offset  constants  are  given  in  Table  A. 2 of  Appendix  A.l.  Many  of  the 
components  are  deemed  non-significant  on  a turning  center  (since  the  turning  center  is  a 
2-axis  machine)  so  that  only  nine  errors  are  included  in  the  final  model  for  the  turning 
center  to  yield  the  error  equations 

Ex  =[£v(z,T)  + £y(xS)]*Zi+£y(zJ)*Z2+[SxXzJ)  + ap(T)*Az]  + Sx(x,T) 

(2) 

Ez=-[£yU,T)  + £y(x,T)]*X4-£y(zJ)*(X2+x)  +Sz(z,T ) 

+ [S.  \x,T)  + ao(T)*  Ax]  - Sz  (s, T) 

Appendix  A. 2 lists  the  program  used  to  symbolically  generate  equations  (1)  and  (2)  using 
the  symbol  manipulator  program  MACS  YMA.  Note  that  the  final  error  compensation 
model  for  the  turning  center  requires  only  two  coordinate  system  offsets.  These  offsets 
are  given  by 


(Z4,X4)  - Tool  offset  from  the  turret. 

(Z2,X2)  - Cross  slide  offset  from  the  carriage  axis. 
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The  notation  T represents  the  vector  of  all  thermocouple  temperatures  given  by 

T = (T0,-Ji9). 


The  following  error  components  on  the  turning  center  are  the  characterization  equations 
and  the  coefficients  used  in  the  current  study  are  given  in  Appendix  B.l.  The  rotational 
displacement  of  all  points  on  the  turning  center  have  a yaw  component  that  is  a function 
of  temperature 

Mach  _ Yaw  = A10  , ((ri0  +T21)1 2)  + A10  2 (3) 

X yaw  is 

ey(x,T)  = Au  x2  + Anx  + (A13  + Mach  _Yaw)  (4) 

Z yaw  is 

ey(z,T)  = A2lz 2 + A22z  + (A23  + Mach  _Yaw)  (5) 

The  X straightness  error  with  Z motion  is 

Sx  \z,T ) = A3iz~  + A22Z  + (A33r39  -i-  A34 ) (6) 

The  X-axis  displacement  is 

3X(X,T)  = A41j:"  + ( A42  + A43T35 )x  + (AUT25  + A45 ) (7) 

The  Z displacement  is 

S_(z,T)  = A^'Z2  + A52z  + (A5}T]2+ Ast)  (8) 

The  Z straightness  of  the  X axis  is 

S.  \x,T)  = A61x‘  + A62x  + (A63T30  + A^)  (9) 

Spindle  growth  is 

S:(s,T)  = AqJn  + A 12T27  + A73T190  + A74  (10) 

The  temperature  of  thermocouple  19  is  taken  at  the  initial  temperature  before  the  machine 
starts  to  warm  up. 

The  squareness  of  the  X-axis  with  respect  to  the  spindle  axis  is 


10 


a0(T)  = \}T0  + (11) 


The  parallelism  error  is 

cc p ( T ) = Aj,Ti9  0 + A92Tl9  0 + Ag3T3]  + 


(12) 
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2.0  Forming  the  Perturbation  Equations 

In  order  to  adjust  the  coefficients  of  (2)  we  form  a set  of  perturbed  equations.  These 
equations  are  built  up  from  perturbing  the  coefficients  of  the  component  error  equations. 
These  are 


£ y (x,T ) + A(fv  (jr,  T))  — | (1  + AAj  i ~ + Ap  (1  + AA)2  )x  + AI3  (1  + AA]3 ) 

+ Al0 , (1  + AAI0 , )((TI0  + T21 ) / 2)  + AI(K2  (1  + AA|0  2 ) 

£ v(z,T)  + A(£  y(z,T))  = A2i  (1  + AA2i  )z~  + A22  (1  + A A22  )z  + A23  (1  + A A23 ) 

+ Ao.i (i + AAio.i  X(P|0  +Tzl)/2)  + A10  2(1  + AAI0  2 ) 

Sx \zj)  + A (Sx \zJ ))  = A3i  (1  + AA31  )z2  + A32(l  + A An)z  + A33(l  + AA33  )T39 
+ A34  (1  + A A34 ) 

Sx  (x,T)  + A (Sx (x,T))  = A41  (1  + AA4I  )x2  + A42  (1  + AA42  )x  + A43  (1  + AA43  )T,5x 
+ A44  ( 1 + A Au  )T35  + A45  (1  + AA45 ) 

3:(z,T)  + A (S.  (z,T))  = A^,  (1  + AA^|  )z~  + A2O  + AA52  )c  + A53(l  + AA^  )7j2 
+ A^l  + AA^4 ) 

S.  \x,T)  + A (S:  \x,T))  = A6,  (1  + AA61  )x2  + A62(l  + AA62  )x  + A63(l  + A \3)T30 
+ A64(l  + AA^) 

Sz  ( s,T ) + A (S.  (s,T))  — A71  (1  + AA7I  )TZ1  + A72(l  + A A72  )T21  + A73(l  + A A73  )7j9  0 

+ A74  (1  + A A74 ) 

ao(T)  + A(ao(T))  = Ag,  (1  + AAg,  )T0  + Ag2(l  + AAg2) 
cq,  CO  + A(tzp  (O)  = A91  (1  + AA^,  )7j9  0 + A92  (1  + AA92  )Ti9  0 + A93  (1  + AA93  )T31 
+ A4  (1 + ) 


Note  that  the  perturbations  are  structured  in  such  a way  that  the  characterization 
coefficients  are  carried  along  into  the  difference  terms.  This  was  done  in  order  to  carry 
the  influence  of  the  characterization  coefficients  along  during  the  fitting  process.  Now, 
to  obtain  the  perturbation  terms  subtract  equations  (4)  through  (12)  from  (15)  to  get 
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A (£y(x,T))  = (Ama2)AAu  + (Ai2a)AA12  + (A13)AA13 

+ (A, 0.1  m*  + T2\  )/2)AA10_,  + (A10  2 )AA10  2 
A(£y(z,T))  = (A2lz2)AA2]  +(A22z)AA22  + (A23)AA23 
+ (A o,i  X(^0  + ^2i  )/2)AA10 , + (A10  2 )AA10  2 
A(J/(z,D)  = (A3lz2)AA31  +(A32z)AA32+(A33r39)AA33 
+ (A34)AA34 

A(t)t(jf,r))  = (A41a2)AA41  + (A42a)AA42  + (A43T35a)AA43 
+ ( a^t^  )aa44  + ( a45  )AA45 

A($;  ( zJ ))  = (A,,  z2  )M5l  + ( A52z)AA52  + (A53ri2)AA53 

+ ( A54  )AAj4 

A (S.  '(x,T))  = (A61a2)AA6i  + (A62a)AA62  + (A637’30)AA63 
+ (A64)AA64 

A (S;(s,T))  = (A7ir22)AA71  +(A72r27)AA72  +(A73r]90)AA73 
+ (A74  )AA74 

A(«o(r))  = (A817’0)AA81  +(A82)AA82 
A(a,(r))  = (^lrl;,0)M„  +(A,,7'I90)AA)2  +(A,37'„)AA„  (16> 

■+  (A94)AA9j 


The  perturbation  form  of  the  kinematic  model  becomes 


A(  Ex ) = [A(£y  (z))  + A(£y  (x)))  * Z4  + A(£y  (z))  * Z2 
+ [A($,  Xz))  + A(ap ) * Az]  + A(«J,  (a)) 


A(£. ) = -[A(£ y (z))  + A(f  v (A))] * X4  - A(£v  (z))  * (X2  + A) 

+ A (S.  (z))  + [A (S:  ’(A))  + A(a0  )*Ax]~  A (Sz  (s)) 

(17) 

In  order  to  relate  the  perturbed  residual  position  errors  in  (17)  to  the  perturbed  error 
component  coefficients  substitute  (16)  into  (17).  Now,  the  problem  becomes  one  of 
determining  the  perturbation  coefficients  from  the  observed  (measured)  residual  error 
vector  (A(EX ),  A (Ey )).  To  do  this  we  create  a system  of  equations  of  the  following  form 
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A (EX)  = CU  (x, z)d} , + C12  (x, z)dn  + C13 (x,  z)d l3 
4~  Cii  (X,  Z ~)d 21  4"  C 22  (X,  Z )d 22  ”1”  ^23  ’ ^)d 23 

+ C31  (x,  z)rf31  + C32  (x,z)d32  + C33  (x,  z)d33  + C34  (x,  z)d34 
+ C4,  (x,  z)d41  + C42  (x,  z)d42  + C43  (x,  zW43  + (x,  z)du  + C45  (x,  z)d45 

+ C91  (x,  z)d9l  + C92  (x,  z)d92  + C93  (x,  z)d93  + C94  (x,  z)^94 
C,0J  (x,  z)d jo  | 4"  ^io,2  2 

A (£’.)  = Du(x,z)dn  + Dl2(x,z)d]2  + Dl3(x,z)dl3 
+ D2I  (x,  z)^21  + £>22  (■*»  Z)^22  + ^>23  (*>  ^23 

+ D51  (x,  z)^51  + £>52  (x,  z)^52  + £>53  (x,  z)rf53  + £>54  (x,  z)d54 
4*  £>61  (x,  z)d b 1 4-  D62 (x,  z)d ^ 2 4-  D63 (x,  z)d(}3  4-  £>(>4  (x,  v, 

+ £>7i  (x,z)rf71  + £>72  (x,  z)d~!2  + D73  (x,  z)d73  4-  D74  (x,  z)d14 
£^81  ^)^81  £>82  ^ )^82 

£>io.l  ZM,o,|  £>10,2  ^)^I0,2 


where  the  Cy  (x,z),£>„  (x,  z),d,;  terms  will  be  defined  below. 

When  (16)  is  substituted  into  (17)  the  unknown  coefficients  are  A Aij . From  the 

numerical  point  of  view  it  will  be  necessary  later  to  scale  these  unknowns  in  equation 
(16)  in  order  to  control  the  magnitude  of  the  elements  in  the  resulting  matrix  that  arises  in 
the  fitting  process.  The  scaled  unknowns  are  the  d coefficients  in  (18).  In  order  to 

construct  them  we  will  introduce  a scaling  array  s(i),  i - l,---,35  and  redefine  the 
perturbation  terms  so  that 
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AA„  = 5(l)*du 
AAI2  =s(2)*d]2 
AA,3  = 5(3)*  dl3 
AA21  =s(4)*d2] 
AA22  =s(5)*d22 
A A23  = s(6)*d23 
AA3I  = s(l)*d2] 
AA32  = 5(8)  * d 22 
M33  =5(9)  *^33 
AA34  -s{\0)*dM 
AA41  = 5(11)*J41 
AA42  = 5(12)  * dx2 
AA43=5(13)*rf43 
= 5(14)*^ 
AA45  = 5(15)  * d45 
AA5I  = 5(16)  * d5i 
AA52  =s(ll)*d52 
AA53  = 5(18)  * d 53 
AA34  = 5(19)*  J54 
AA61  =5(20)*  d6i 
AA62  = s(2l)*  d62 
AA63  =5(22)*  db3 
AAm  = 5(23)  * dM 
AA7I  = 5(24)  *d7l 
AA12  - s(25)*  d12 
AA„  = s(26)*  d13 
AA74  = 5(27)  * d74 
AAgl  = 5(28)  *dgl 
AAg2  =s(29)*dS2 
AA,,  = 5(30)  *d9i 

AAg2  = 5(31)  *d92 

AA^3  = 5(32)  * d93 
A A^  =5(33)*c(94 
AAio.i  = 5(34)*  JI01 
AA|o  2 =5(35)*  dl0  2 


(19) 
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where  d tJ  are  new  perturbation  coefficients.  The  scaling  is  necessary  in  order  to 

guarantee  that  terms  in  the  fitting  process  matrix  are  all  approximately  the  same  order  of 
magnitude  and  are  not  too  large.  The  scaling  array  is  determined  in  such  a way  that  the 
maximum  value  in  each  column  of  the  matnx  J is  one.  This  is  done  for  each  column  by 
first  finding  the  maximum  of  the  absolute  values  of  the  elements  in  each  column.  If  the 
maximum  is  zero  then  the  scale  for  that  column  is  set  to  one.  Otherwise  the  scale  is  set  to 
the  reciprocal  of  the  maximum.  For  a discussion  of  scaling  and  its  application  to 
problems  of  inversion  see  [25]. 

Returning  to  the  substitution  of  (16)  into  (17)  we  further  introduce  (18)  into  (17)  and 
rearrange  terms.  Now,  define  the  following  coefficients 


(A,Z) 

= 

(A, 

* 5(1)  * Z4  )x2 

c 
^ 1 2 

(x,Z) 

= 

(A,! 

* 5(2)  * Z4  )x 

C13 

(x,z) 

= 

^13 

*s(3  )*Z4 

C2I 

(a,z) 

= 

(A, 

* 5(4)  * (Z4  4 

■z3 

c 22 

(x,z) 

= 

(A22  * 5(5)  * (Z4  +Z2)). 

c2i 

(x,z) 

= 

^23 

* 5(6)  * (Z4  + 

z2) 

C3i 

(a,z) 

= 

^31 

*5(7)  * Z2 

Ci2 

(a,z) 

= 

^32 

* 5(8) * z 

Ci3 

(a,Z> 

= 

^33 

*T39  *5(9) 

c34 

(a,z) 

= 

^34 

* 5(10) 

C4I 

(a,Z) 

= 

*4, 

* 5(1  1)  * x 2 

C42 

(x,z) 

= 

^42 

* 5(12)  * Jt 

C43 

(a,z) 

= 

^43 

*r35  *5(13)* 

X 

cu 

(x,z) 

= 

^44 

*r35  *5(14) 

(x,z) 

= 

A45 

*5(15) 

C9I 

(A,  z) 

= 

A9| 

*rI20  * 5(30) 

* Az 

c 

^ 92 

(a  ,Z) 

= 

A92 

*T19.0*5(31) 

* Az 

c 

^93 

(a,  z) 

= 

A93 

*r3l  *5(32)* 

Az 

r 

94 

(a  ,z) 

= 

A94 

* 5(33) * Az 

C10,  (a,  z)  = A10a  * (2 Z4  + Z2  )((T10  + r2I ) / 2)  * s( 34) 
c, 0.2  (a,  z)  = Al0  2 * (2Z4  + Z2 ) * 5(35) 


and 
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Du(x,z)  = -(Au  *s(l)*XA)x2 

Dn  (x,  z ) = -(A,  2 * 5(2)  * X4)x 

Dl3(x,  z)  = -A13  * 5(3)  * X4 

D2i(jc,z)  = -(A21  *s(4)*(X,  +X2  + x))z2 

D22  (x,z)  = — (A22  * 5(5)  * (X4  + X2  + x))z 

D23(x,z)  = -A23  * 5(6) * (X4  + X2  + x) 

D5](x,z)  = A5I  *5(16)*z2 
D52(x , z)  = A52  * 5(17)*  z 
D52(x,  z)  = A53  *Ti2  *5(18) 

D^ix,  z)  = Am  *5(19) 

D6l(x,z)  = A6i  *5(20)*  a2 
D62(x,z)  = A62  *5(21)*  x 
D63(x,z)  = A63  *r30  *5(22) 

Dm(x,z)  = Am  *5(23) 

D7l(x,z)  = -A7l*T22  *5(24) 

D72(x,z)  = -A72  *T27  * 5(25) 

D73  (x,  z)  = — A73  * r]9  0 * 5(26) 

D7i(x,z)  = -A74  * 5(27) 

Dgl  (a, z) = A81  * T0  * 5(28)  * Aa 
D82(a,  z)  = A82  * 5(29)*  Ax 

dio.iUz)  = -A,o.i  *(X2  +2X4  + a)*((Ti0  +r21)/2)*5(34) 
7)i0i2(^^)  = -Ai02  *(X2  +2X4  + a) * 5(35) 


By  substituting  equations  (20)  and  (21)  into  (17)  we  get  (18),  which  is  the  system  we 
sought  to  form.  System  (18)  itself  cannot  be  directly  solved  for  the  unknowns,  but  if  we 
introduce  the  residual  errors  measured  at  a set  of  spatial  points  on  the  part  then  it  will  be 
feasible  to  solve  for  the  unknowns  in  a least  squares  manner.  In  the  next  section  we  form 
the  system  of  equations  needed  to  solve  for  the  unknowns.  There  are  35  unknowns,  each 
associated  with  a column  of  the  system  matrix.  However,  not  all  of  the  coefficients 
appear  in  both  (20)  and  (21)  simultaneously.  These  missing  coefficients  lead  to  blocks  of 
zeroes  in  the  matrix  defined  in  the  next  section. 
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3,0  Forming  and  Solving  the  System  Matrix 

We  form  the  system  of  equations  to  solve  for  the  unknown  dt]  's  in  equation  (18)  as 
follows.  Let  there  be  n points  on  the  part  at  which  residual  errors  are  to  be  computed 
from  CMM  measurements.  Denote  the  points  by  (jq , z, ) for  i = 1,  n.  Then  set 


J,.\  = 
1.2  = 

J 


i,3 


J.,  = 


1.4 


I. 5  = 

J,  6 = 

J,,= 

J i.  8 — 

J,  9 = 

A. 10 

J,n  = 

^ i.l2 

= 


1.14 


J 

1.15 

I, 

J 


i.30 


J 


i,31 


A. 32  ~ 
•^1,33  = 


i.34 


i,35 


C\ 1 (-*o+]-i  ’ *“o+I-i  ) 
^12  (**0+1 -i  ’ ^o+l-i  ) 
^13  (-*o+l-i  » ^o+l-i  ) 
^21  (Xn+\-i  ’ *“o+l-i  ) 
^22  (^o+l-i  ’ *L+l-i  ) 
^23  (*o+l-i  ’ ^o  + l-i  ) 

— C3I  (-Ti+l-i  , ^n+i-,  ) 

= C32  (-^o+l-i  ’ ^o+l-i  ) 

— C 33  (■*„+]_,  , £„+]_,  ) 

— (^34  (-*„+!-,  i ^n+l-i  ) 
^41  (**o  + l-i  ’ *“0+1-1  ) 

— C 42  (■*„+]_,  , Zn+]_,  ) 

— C43  , Z„+\-i  ) 

— C44  (-^o+l-i  ’ *L+l-i  ) 

— L45  (jtn+|_,- , Zn  + 1_,  ) 

= 0,k  = 16,  ••*,29 

— C91  (-Tj+i-,  ’ *“n+i-i ) 

— ^92  (-^o+i-i  ’ ^o+l-i  ) 

— C93  (jc„+1_1-  , zn+1_J ) 

— C94  (xf;+,_. , zn+I_, ) 

— Cl0  j (-X„+1_, , v!n+|_,  ) 

~ ^10,2  (*o+l-i  ’ ^0+1-1  ) 


(22) 


In  C91  to  C94  the  value  of  A z is  computed  as  Az  = zn+l_,  - z„_,  for  / = !,*••,«  - 1 . For 
i = n,  A z = z2~zr 

For  i = n+l,...,2n  the  matrix  elements  are  written  as 
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0,  1 (-*2,1+1-,  ’ 

2/i+l—j 

) 

1.2 

= 0,2' 

(*2 /,+!-, 

’ ^2/i+l-i 

) 

ly 

0i3(^2n+l-i  ' 

> ^2«+l-, 

) 

J,A 

= O+l1 

(-*2/i+l-i 

’ ^ 2/,+1-z 

) 

J,  5 

= Z)22 

("*2/i+l-i 

’ ^ 2/i+ 1— i 

) 

J,  6 

= F)+3 ' 

(-*2/1+!-: 

’ *”2/i+l-< 

) 

*/,,* 

= 0 ,k 

= 7,--- 

,15 

*/„,6 

= 051 

(-*2n+l-; 

» ^2/i+l- 

,) 

A,n 

= 052 

(•*2/!+,-/ 

’ *”2/i+l- 

,) 

*^,.18 

= 053 

(•*2/1  + 1-/ 

’ ^ 2/i+l- 

,) 

^ i,l9 

= 054 

("*2/1+1-, 

’ *”2/i  + l- 

,) 

J i,  20 

= 061 

(-*2/1  + 1-/ 

’ *”2/i+l- 

,) 

A, 21 

= 062 

(-*2/1+1—, 

’ ^ 2/i+ 1 — 

,) 

22 

= 063 

("*2/1+1-/ 

; * *”2/i+l- 

,) 

*^,.23 

= 064 

("*2/i+l-, 

’ *”2/i+l- 

,) 

A. 24 

= 07, 

(•*2/1  + 1-, 

5 ^ 2/1+1— 

,) 

^ i,25 

= 072 

(-*2/1+1-, 

- ’ *”2/i+l- 

,) 

J i.26 

= 073 

(■*2,1+1-, 

’ ^ 2/i+l  — 

,) 

J i,21 

= 074 

(-*2n+l-i 

! ’ *”2/1  + 1- 

,) 

J i,  28 

= 081 

(-*2/1+1-, 

’ 2/i  + l- 

,) 

^ i,29 

= 082 

(-*2/1  + 1-, 

’ ^ 2/i+l- 

,) 

J.,k  ■ 

= o. 

k=  30 

,--•,33 

Im 

010,1  (-*2/1+1- 

-,  ’ ^2/i+l 

-,) 

J i,35 

— 0,0 

,2  (-*2/1+1 

-,  ’ ^ 2/1+! 

In  Dgl  and  D82  the  value  of  Ax  is  computed  as  Aj:  = jc2n+l_(  - x2n_t  for 
t = « + 1,  • • • ,2/2  - 1 . For  i = 2n,  Ax  = x2  - xr 


Finally,  for  ease  of  notation,  set 

V(i)  = AEnx+'-,i  = !,■■■, n 
V (i)  = AE2n+l~‘  ,i  = n + !,■■■  ,2n 


(24) 


and 
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U(l)  = du 
U (2)  = dn 
U (3)  = J13 
U (4)  = d2x 
U (5)  = d22 
U (6)  = d23 

uv)  = d3l 

U (8)  = di2 
U(9)  = d» 

um  = dM 

U(ll)  = dM 
U(l2)  = d42 
U(l3)  = di3 
U (14)  = du 
U (15)  = d45 

um  = d5] 

U (17)  = d52 

um  = d52 

U (19)  = d54 
U (20)  = d6l 
U (21)  = d62 
U{22)  = d63 
U (23)  = dM 
U(24)  = dlx 
U (25)  = d72 
U (26)  = d12 
U (27)  = d14 
U (28)  = dix 
U(  29)  = di2 
U (30)  = d9l 
U(3l)  = d92 
U (32)  = d93 
U(33)  = d94 
U (34)  = d]OX 
U (35)  = dl02 


(25) 
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We  now  have  the  notation  set  up  to  form  the  system 

JU  = V (26) 

where  J is  a matrix  with  2n  rows,  since  there  are  n probe  points  each  with  an  x and  a z 
error,  by  35  columns,  one  for  each  scaled  perturbation  coefficient.  Note  that  if  any 
residual  error  in  equation  (24)  is  zero  then  the  associated  row  of  J is  set  to  zero. 

If  J were  a square  matrix  and  non  singular  then  the  solution  would  be  simply 

(27) 

but  in  general  J is  an  overdetermined  system  with  more  rows  than  columns.  Furthermore, 
since  many  of  the  columns  might  be  dependent  on  each  other  a standard  least  squares 
approach  will  likely  fail.  For  this  reason  we  apply  a matrix  called  the  generalized  or 
pseudo  inverse  and  compute 


U = J + V.  (28) 

The  details  involved  with  computing  the  pseudo  inverse  are  given  in  Appendix  E and  in 
[24]. 

Once  U has  been  computed  the  A A coefficients  can  be  determined  from  (19)  and  the 
updated  coefficients,  given  by 

A/l  + AA,,).  (29) 
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4.0  Application 


In  this  section  an  application  of  the  inversion  process  described  in  the  previous  sections 
will  be  given.  The  part  shown  in  Figure  2 was  machined.  It  consisted  of  two  concentric 
cylinders.  The  smaller  cylinder  has  a nominal  diameter  of  50.8  mm  (2  in),  and  the  larger 
one  has  a nominal  diameter  of  152.4  mm  (6  in).  To  compute  the  process  errors,  the  part 
was  probed  with  a touch  trigger  probe  at  35  points  while  still  secured  on  the  machine 
tool.  For  a discussion  of  on-machine  probing  see  Bandy  and  Gilsinn  [17]. 


Figure  2:  Two  Cylinder  Turned  Part 


In  order  to  use  on-machine  touch  trigger  probe  readings  in  determining  machine  tool 
errors,  various  coordinate  systems  have  to  be  linked.  For  example,  the  turning  center 
error  models  are  based  on  calibrated  error  component  equations  that  describe  error 
motions  at  points  in  the  reference  coordinate  system  of  the  turning  center.  When  parts 
are  turned,  however,  points  on  the  part  are  identified  relative  to  a part  coordinate  system. 
Therefore  the  part  coordinate  zero  needs  to  be  located  relative  to  the  turning  center 
reference  coordinate  system.  In  Figure  2 the  part  coordinate  zero  is  set  as  the  center  point 
of  the  top  face  of  the  small  cylinder.  The  positive  Z-axis  is  taken  vertically  and  the 
positive  X-axis  is  taken  to  the  right.  In  order  to  link  the  part  zero  point  to  the  turning 
center  reference  coordinate  system  two  sets  of  data  are  read.  First  the  probe  offsets  in  the 
X and  Z directions  are  obtained  from  a database  file.  These  offsets  are  distances  from  a 
point  on  the  turning  center  tool  turret  called  the  tool  setting  point.  As  the  turret  moves 
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the  location  of  the  tool  setting  point  with  respect  to  the  reference  coordinate  system  can 
be  queried  so  that  when  the  touch  trigger  probe  is  brought  to  the  part  zero  point  the 
location  of  the  part  zero  relative  to  the  reference  coordinate  system  can  be  computed.  As 
a result  any  part  coordinate  can  be  identified  with  a reference  coordinate.  Therefore, 
when  points  on  a part  are  identified  by  their  part  coordinates,  error  estimates  at  these 
points  can  then  be  computed  by  identifying  the  part  coordinates  with  their  machine 
reference  coordinates.  For  a further  discussion  of  computing  machine  offsets  see 
Appendix  C. 

A coordinate  measuring  machine  (CMM)  is  used  to  measure  the  part  in  Figure  2 at  the 
same  part  coordinates  that  were  probed  on  the  turning  center.  This  is  done  by 
establishing  a working  reference  plane  on  the  CMM  along  the  top  surface  of  the  small 
cylinder  and  setting  the  reference  plane  zero  at  the  same  point  as  the  part  zero  is  set  on 
the  turning  center.  Then,  the  35  points  probed  on  the  turning  center  are  probed  at  the 
same  coordinates  on  the  CMM.  Error  measurements  are  made  normal  to  the  surface 
features  probed.  Realizing  that  there  is  statistical  variation  in  the  measurements,  the 
CMM  is  programmed  to  measure  the  35  points  along  the  surface  at  8 angles  around  the 
part  in  Figure  2.  The  error  measurements  at  these  equally  spaced  angles  are  then 
averaged. 


The  estimation  of  the  new  coefficients  for  the  component  error  equations  in  the  turning 
center  error  model  proceeds  by:  1)  reading  the  old  coefficients  from  a database,  2) 
reading  the  thermocouple  temperatures  on  the  turning  center  when  the  part  was 
manufactured,  3)  reading  the  initial  turning  center  temperatures  before  machining,  4) 
reading  the  nominal  points  probed  on  the  part  and  the  averaged  CMM  error 
measurements.  The  turning  center  errors  are  computed  by  the  methods  of  Appendix  D 
and  the  system  matrix  J is  formed  as  in  Section  4.  Next,  the  predicted  turning  center 
errors  are  computed  using  the  old  coefficients  and  the  difference  between  the  CMM 
measured  errors  and  the  predicted  turning  center  errors,  as  in  equation  (D5).  This 
produces  the  residual  errors  at  the  probed  points.  The  pseudoinverse  of  the  system 
matrix,  J + , is  applied  to  the  residual  errors  to  produce  the  scaled  correction  terms  dl} . 

The  coefficient  correction  terms  are  computed  by  rescaling  the  terms  in  (19).  Finally  the 
new  coefficients  are  computed  by  adding  one  to  the  rescaled  correction  terms  and  then 
multiplying  by  the  old  coefficients  as  in  equation  (29). 

The  computed  results  are  shown  below  in  figures  3 and  4.  Of  the  35  points  chosen  on  the 
surface  for  analysis,  eight  were  selected  on  the  top  surface  of  the  small  cylinder  and  nine 
each  were  selected  on  the  small  cylinder  side,  the  large  cylinder  top  and  the  large  cylinder 
side.  Figure  3 shows  the  results  of  plotting  the  averaged  CMM  error  measurements  in  the 
X direction  along  with  the  predicted  machine  errors  using  the  updated  coefficients.  The 
solid  lines  are  the  predicted  errors  and  the  dashed  lines  are  the  averaged  CMM  error 
measurements  at  the  points  along  the  small  cylinder  (points  9 to  17)  and  along  the  large 
cylinder  (points  27  to  35).  Figure  4 shows  a similar  plot,  but  for  the  errors  in  the  Z 
direction.  The  errors  along  the  top  of  the  small  cylinder  are  plotted  for  points  1 through  8 
and  those  along  the  top  of  the  large  cylinder  are  for  points  18  through  26.  The  maximum. 
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mean  and  the  standard  deviation  of  the  differences  between  the  measured  and  predicted 
errors  are  given  in  Table  1.  The  new  coefficients  are  given  in  Appendix  B2. 


Error 

Difference  (mm) 

Maximum 

Mean 

Standard  Deviation 

X Direction 

0.0041 

0.0014 

0.0010 

Z Direction 

0.0025 

0.0008 

0.0007 

Table  1:  This  shows  the  Maximum,  Mean  and  Standard  Deviations  of  the 
differences  in  the  errors  plotted  in  figures  3 and  4. 


X Direction 


Point  Number 


Figure  3:  Comparing  the  Measured  Turning  Center  Errors  (dashes)  with  the 
Predicted  Errors  (solid)  based  on  Updated  Turning  Center  Model  in  the  X 

Direction. 
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Z Direction 


Figure  4:  Comparing  the  Measured  Turning  Center  Errors  (dashes)  with  the 
Predicted  Errors  (solid)  based  on  Updated  Turning  Center  Model  in  the  Z 

Direction. 
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5.0  Discussion 


Machine  tool  error  models  are  often  developed  and  calibrated  by  using  such  techniques 
as  laser  interferometry  and  capacitance  gauges  to  measure  component  errors  such  as  scale 
and  straightness  errors.  Since  a machine  tool’s  various  parts  such  as  bearings  cannot  in 
general  be  warmed  up  to  functional  conditions  experienced  during  machining,  another 
method  must  be  developed  to  modify  the  calibrated  machine  error  models  to  represent 
errors  experienced  under  cutting  conditions.  This  report  proposes  a method  that  can  be 
used  to  modify  the  coefficients  of  a calibrated  model  so  that  it  can  be  used  during 
machining  operations  to  compensate  for  machine  tool  errors. 

The  method  employed  here  is  a modified  least  squares  technique  that  takes  into  account 
the  fact  that  the  normal  equations  generated  during  the  least  squares  process  may  not  be 
solvable  in  the  usual  sense.  For  that  reason  a generalized  inverse  is  applied  to  the  normal 
equations  to  generate  the  best  solution  in  the  sense  of  Appendix  E.  This  is  a well  known 
tool  for  solving  least  squares  problems  and  a thorough  discussion  can  be  found  in  [24], 

The  results  shown  in  Table  1 and  Figures  3 and  4 clearly  show  that  the  method  produces 
an  adjusted  error  model  that  accurately  reproduces  the  measured  data  used  to  update  the 
error  model  for  the  turning  center  used  during  the  study.  Appendices  B1  and  B2  show 
that  the  method  significantly  changes  the  coefficients  to  the  basic  form  of  the  model.  A 
reader  should  recognize  that  the  model  fits  are  adequate  for  the  range  of  the  data.  They  do 
not  guarantee  that  the  updated  model  could  be  used  to  compensate  errors  on  turned  parts 
dissimilar  to  the  one  used  in  the  study.  Further  research  must  address  the  question  of  how 
to  extend  the  model  adjustment  procedure  described  in  this  report  to  hold  for  families  of 
similar  parts  or  for  other  classes  of  parts.  Further  studies  will  also  emphasize 
manufacturing  parts  under  varying  temperature  conditions.  Parts  will  also  be  designed  to 
test  different  features  placed  throughout  the  workspace 

Since  the  intent  of  the  machine  tool  error  model  is  to  be  used  for  real  time  error 
compensation  no  attempt  is  made  to  fit  model  forms  beyond  cubic  spatial  terms.  This  is 
because  past  experience  has  shown  that  cubic  polynomials  can  be  evaluated  and  their 
results  passed  to  a machine  tool  to  correct  tool  position  within  the  real  time  constraints  of 
the  cutting  process.  Furthermore,  higher  order  polynomials  can  introduce  unwanted 
oscillations  between  fitting  points. 

Again  the  authors  emphasize  that  this  is  not  a report  that  claims  that  using  singular  value 
decomposition  is  the  only  tool  for  modifying  or  developing  a machine  tool  error 
compensation  model.  However,  given  the  initial  error  component  model  forms  on  which 
to  base  the  turning  center  error  compensation  model,  the  authors  felt  that  singular  value 
decomposition  was  a sufficiently  robust  tool  to  use  for  the  current  study. 
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APPENDIX  A 


A.l  Developing  the  Kinematic  Model  for  a Turning  Center 

In  this  appendix  the  homogeneous  transformations  for  the  coordinate  systems  modeling 
the  turning  center  components  along  with  the  final  error  model  will  be  given.  For  a 
detailed  discussion  of  the  theoretical  development  of  the  transformations  see  [3-5],  The 
turning  center  that  is  modeled  in  this  section  is  shown  in  Section  2,  Figure  1.  Table  A1 
lists  thermocouples  along  with  their  position  descriptors. 

The  first  step  in  the  model  development  is  the  establishment  of  a reference  coordinate 
system.  On  the  turning  center  this  will  be  taken  as  the  or  absolute  coordinate  system,  the 
origin  of  which  is  located  at  the  farthest  travel  to  the  left  along  the  z-axis  and  the  farthest 
travel  vertically  alone  the  x-axis.  For  the  turning  center,  modeled  in  this  study,  the 
absolute  coordinate  system  origin  aligns  itself  with  a point  that  is  located  at  the  upper 
right  rear  of  the  tool  slot  as  one  faces  the  turret,  when  the  turret  is  brought  to  Home 
position.  This  point  is  also  called  the  tool  setting  point. 

Let  (Z, , X, ) be  the  offset  of  the  origin  of  the  z-slide  coordinate  system  from  the 
reference  coordinate  system.  Upper  case  letters  are  constants  representing  coordinate 
system  offsets  from  each  other,  whereas  lower  case  axis  letters  represent  variable  values. 
Thus  the  origin  of  the  z-slide  axis  system  is  set  Z\  units  along  the  z-axis  from  the  origin 
of  the  reference  coordinate  system  and  X\  units  along  the  x-axis  from  the  origin  of  the 
reference  coordinate  system.  The  motion  of  the  z-slide  coordinate  system  (also  called  the 
carriage  coordinate  system)  relative  to  the  machine  coordinate  system  is  described  by  the 
homogeneous  transformation 

' 1 -€Z(Z,T)  £y(Z,T) 

°T  = 1 ~£AzJ) 

— £ y (z,T)  £x(z,T)  1 

0 0 0 

£_(z,T)  is  the  roll  error  of  the  carriage  given  z motion  as  a function  of  the 
temperature  state,  T,  of  the  machine. 

£x(z,T ) is  the  pitch  error  of  the  carnage  given  z motion  and  temperature 
state  T. 

e ( z,T ) is  the  yaw  error  of  the  carriage  given  z motion  and  temperature 

state  T. 

Sx(z,T)  is  the  x-displacement  error  given  z motion  and  temperature  state  T. 

Sy  ( z,T ) is  the  y-displacement  error  given  z motion  and  temperature  state  T. 

8 , (z,T)  is  the  z scale  error  given  z motion  and  temperature  state  T. 


Xl+8x(z,T) 

Sy(z,T) 

Zx  + z + 8z(z,T) 

1 


(A.l) 
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Sx(z,T)  = SxXz,T)  + ap(T)*Az 


(A. 2) 


where 

8x\zJ) 

ap(T) 


Az 

T 


is  the  x-straightness  error  of  the  carriage  given  z motion  and 
temperature  state  T. 

is  the  parallelism  error  between  the  z-motion  and  axis  average  line 

of  the  spindle  given  temperature  state  T. 
is  the  incremental  z motion. 

temperature  state  of  the  machine  and  is  taken  to  be  the  vector 
T = (T0,---,r39) , of  thermocouple  temperatures. 


Let  (Z2,X2)  be  the  offset  of  the  origin  of  the  x-slide  axis  system  (also  called  the  cross 
slide)  from  the  origin  of  the  z-slide  axis  system  then  the  transformation  representing  the 
motion  of  the  x-slide  motion  relative  to  the  z-slide  is  given  by 


ZT  = 


1 

£;(x,T) 

-ey(x,T) 

0 


-e:(x,T) 

1 

£x(x,T) 

0 


£y(xJ)  X2+x  + Sx(x,T)^ 
-£x(x,T)  Sy(x,T) 

1 Z2  +S.(x,T) 

0 1 


(A. 3) 
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Thermocouple  Locations 

Number 

Location 

0 

Bottom  of  X-Scale  (Glass) 

1 

Top  of  X-Scale 

2 

Coolant  Tank 

3 

Not  Used 

4 

Not  Used 

5 

Not  Used 

6 

Not  Used 

7 

Right  of  Z-Scale 

8 

Right  Center  of  Z-Scale 

9 

Left  Center  of  Z-Scale 

10 

Left  of  Z-Scale 

11 

Top  of  X-Way 

12 

Bottom  of  X-Way 

13 

Top  of  X-Head 

14 

Bottom  of  X-Head 

15 

Bottom  Left  of  Z-Slide 

16 

Top  Left  of  Z-Slide 

17 

Bottom  Right  of  Z-Slide 

18 

Top  Right  of  Z-Slide 

19 

Hydraulic  Tank 

20 

Left  End  of  Lower  Z-Way 

21 

Left  End  of  Upper  Z-Way 

22 

Right  End  of  Lower  Z-Way 

23 

Right  End  of  Upper  Z-Way 

24 

Lower  Front  of  Spindle  Head 

25 

Lower  Rear  of  Spindle  Head 

26 

Upper  Front  of  Spindle  Head 

27 

Upper  Rear  of  Spindle  Head 

28 

Left  of  Top  of  Spindle  Head 

29 

Middle  of  Top  of  Spindle  Head 

30 

Right  of  Top  of  Spindle  Head 

31 

Bottom  Left  of  Bed 

32 

Top  Left  of  Bed 

33 

Bottom  Right  of  Bed 

34 

Top  Right  of  Bed 

35 

Near  X-Drive  Motor  Shaft  Bearing 

36 

Left  Z-Ballscrew  Bearing 

37 

Right  Z-Ballscrew 

38 

X-Ballscrew  Housing 

39 

Z-Ballscrew  Nut 

Table  Al:  Thermocouple  Locations  on  the  Turning  Center. 
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where 


e,(x,T) 

is  the  roll  error  of  the  cross  slide  with  x motion  at  temperature  state 

T. 

e.(x, T) 

is  the  pitch  error  of  the  cross  slide  with  x motion  at  temperature 
state  T. 

Cy{X,T) 

is  the  yaw  error  of  the  cross  slide  with  x motion  at  temperature 

state  T. 

S,(x,T) 

is  the  y displacement  of  the  cross  slide  with  x motion  at 
temperature  state  T. 

S,(x,T) 

is  the  x scale  error  of  the  cross  slide  with  x motion  at  temperature 

state  T. 

S.( x,T) 

is  the  z displacement  error  of  the  cross  slide  with  x motion  at 
temperature  state  T. 

S:(xJ)  = S:XxJ)  + ao(T)*Ax  (A.4) 

where 

Sz\x,T) 

is  the  z straightness  of  the  cross  slide  with  x motion  at  temperature 
state  T. 

&JT) 

is  the  orthogonality  error  between  the  x motion  and  the  axis 
average  line  of  the  spindle  at  temperature  state  T. 

Ax 

is  the  incremental  x motion. 

Let  the  turret  and  x-slide  coordinates  overlay  so  that 


' 1 

£z(tJ) 
~ £y(t,T ) 


-et(t,T)  ey(t,T ) S,(t,T)' 

1 -c,(t,T)  5y(t,T) 

e,(t.T ) 1 SAt,T) 

0 0 1 


where,  due  to  misalignments  of  the  turret  on  its  rotational  axis, 

£(t,T),£(t,T),£.(t,T ) are  the  rotational  errors  in  station-to-station  motion  of 
the  tool  turret  at  temperature  state  T. 

Sr(t,T),S(t,T),S.(t,T ) are  the  translation  errors  in  station-to-station  motion  of 

x y c 

the  tool  turret  at  temperature  state  T. 
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A turret  on  the  lathe  has  the  capability  of  holding,  for  example,  various  cutters,  tools  and 
on-machine  measurement  devices.  Station-to-station  motion  here  refers  to  rotation  of  the 
turret  from  one  cutter  or  tool  position  on  the  turret  to  another. 


Let  (Z4 , X4 ) be  the  offset  of  the  cutter  coordinate  system  from  the  turret  coordinate 
system.  Since  the  cutting  tool  coordinate  frame  is  assigned  to  the  tool  tip  (a  point),  the 
rotation  errors  have  no  effect.  The  errors  involved  are  due  in  part  to  the  fact  that  an 
actual  tool  is  in  error  from  its  nominal  dimension  and  in  part  to  thermal  effects  during  the 
cutting  operations.  The  cutting  tool  error  matrix  relative  to  the  turret  is  given  by 


^ 1 0 0 X4+£,(c,7T 

0 1 0 Sy(c,T) 

0 0 1 ZA+S:(cJ) 

0 0 0 1 

V J 


(A. 6) 


where 

Sf(c,T),Sv(c,T),SAc,T),  are  the  incremental  tool  dimension  changes  in  the  x,  y 
and  z directions  at  temperature  state  T. 


Let  (Z5,  X5 ) be  the  offset  of  the  spindle  coordinate  system  from  the  machine  coordinate 
system.  In  the  generic  turning  center  model,  the  spindle  has  unrestricted  rotation  about 
the  z axis  so  that  £.  (s,T)  - 0 . The  spindle  errors  are  due  to  either  axis  of  rotation  errors 
of  the  spindle  or  errors  resulting  from  misalignments  of  the  chuck.  The  spindle  error 
matrix  relative  to  the  machine  coordinate  system  is  given  by 


( 


V 


1 

0 


£y(s,T) 

0 


0 

1 


£x(s,T) 

0 


ey(s,T)  X5+SxU,T)'' 
~£As>T)  Sr(s,T ) 

1 Zt+SAs,T ) 


(A.7) 


where 

£x(s,T) 
£y(sJ) 
Sx(s,T) 
Sy(s,T ) 
Sz(s,T) 


is  the  tilt  of  the  spindle  in  the  non-sensitive  direction,  i.e.  yaw 
motion  about  the  x axis,  at  temperature  state  T. 
is  the  tilt  of  the  spindle  in  the  sensitive  direction  or  a pitch  motion 
about  the  y axis  at  temperature  state  T. 

is  the  translations  of  the  spindle  in  the  sensitive  or  x direction  at 
temperature  state  T. 

is  the  translations  of  the  spindle  in  the  non-sensitive  or  y direction 
at  temperature  state  T. 

is  the  z-axial  motion  of  the  spindle  at  temperature  state  T. 


33 


Let  (z(w),  jt(w))  be  the  desired  point  on  the  work  piece  to  which  the  tool  is  to  be  brought. 
Finally,  for  the  workpiece  assume  no  rotational  errors.  The  workpiece  error 
transformation  is  given  by: 


"10  0 x(w)  + Sx(w,T)^ 
0 1 0 Sy(w,T) 

0 0 1 z(w)  + Sz  (w,T) 

0 0 0 1 


(A. 8) 


where 

Sx(w,T),S  y(w,T),S.(w,T)  are  the  changes  of  the  coordinates  of  the  point  on  the 

workpiece  to  which  the  tool  is  to  be  applied  in  the  x,  y and  z directions  at 
temperature  state  T. 


To  determine  (z(w),jc(w))  at  the  work  point  one  equates  the  ideal  transformation  (i.e. 
without  errors)  products  of  the  two  chains 


or  , ° (A>9) 


and  get 


a:(w)  — X,  + X + X4  + x — X5 
z(w)  = Z]+Z2+Zi+z-Z5 


(A. 10) 


The  offsets  are  given  in  Table  A. 2 for  the  turret  lathe.  The  coordinate  systems  are  not 
placed  arbitrarily  but  are  placed  on  the  machine  tool  components  in  such  a way  that  they 
are  related  to  the  measurement  characteristics  of  the  metrology  devices  used  to  calibrate 
the  error  components. 


X 

Z 

1 

-148 

-43 

2 

30 

-26 

3 

0 

0 

4 

-35 

30 

5 

-386 

0 

Table  A.2:  Coordinate  System  Offsets 
in  mm. 
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The  resulting  generic  error  matrix  can  be  computed  as 


E =(sT~l  )°Tz  Tx  % % (ATI) 

A symbolic  manipulation  package,  such  as  MACSYMA,  can  be  used  to  generate  the 
inverse  matrices  and  compute  the  final  product  E (see  Appendix  A. 2.).  The  final 
displacement  errors  of  the  tool  from  the  workpiece  are  given  by: 

Ex=£y(sX*Z5+(ey(zX  + £y(x,T)  + £y(t,T)-£y(s,T))*Z4 
+ (ey(z,T)-£y(s,T))*Z2 -£y(s,T)*(Zi+z)  + Sx(z,T)  + Sx(x,  T) 

-Sx(w,T)  + Sx  (r,  F)  - Sx  (s,T)  + 8X  (c,  T) 

Ey  = -ex  (5,  T)  * Z5  + (~ex  (z,T)~£x(x,T)~  £x  ( t,T)  + £x  (s,  T))*Z 4 

+ (-£x(z,T)  + £x(s,T))*Z2  +ex{s,T)*(Zi  +z)  + (£z(z,T)  + £:(x,T)  + £:(tS))*X4  (A'12) 

+ £z(zX)*(X2+x)  + Sy(z,T)  + Sy(x,T)-Sy(w,T)  + Sy(tX)-Sy(s,T)  + Sy(c,T) 

E:  =-ey(s,T)*Xs+(-£y(z,T)-£y(x,T)-£y(t,T)  + £y(sX))*X4 
+ (-£y(z,T)  + ey(s,T))*(X2  + x)  + £y(s,T)*  X,  +Sz(z,T)  + Sz(x,T) 

- 5.  (• w,T ) + Sz  ( tj ) - 8.  (5,  T)  + 8Z  (cj) 

These  equations  give  the  generic  error  formulas  for  a turning  center  with  its  absolute 
reference  coordinate  system  and  the  tool  setting  point  assumed  aligned  at  “cold  start”  in 
the  “home”  position. 

There  are  only  two  equations  in  the  final  model.  Errors  in  the  y direction  are  not  included 
since  they  are  in  the  non-sensitive  direction.  The  nonsensitive  direction  in  turning  centers 
is  the  direction  perpendicular  to  the  plane  in  which  the  two  machine  slides  move.  The 
only  error  components  included  were  those  found  to  be  significant  and  not  in 
nonsensitive  directions.  Only  seven  errors  are  included  in  the  final  model  for  the  turning 
center  to  yield  the  error  equations 

E,=(£)(z,T)  + ey(x,T))*Z,+ey(z,T)*Z2+(Sx’U,T)+ap(T)*Az)  + S,(x,T) 

(A. 13) 

E.  = -(ey(z,T)  + £y(x,T))*  X4  - ey(z,T)*(X2  + x)  +Sz(z,T) 

+ (SAx,T)  + ct„(T)*Ex)-SAs,T) 

The  MACSYMA  script  used  to  generate  these  equations  is  given  in  Appendix  A. 2. 
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A.2  MACSYMA  Script  for  Kinematic  Model  Generation 

This  appendix  contains  the  listing  of  a command  script  for  the  symbol  manipulator 
language  MACSYMA.  The  script  can  be  loaded  into  MACSYMA  by  using  the  BATCH 
command.  Note  that  the  small  variable  arithmetic  is  taken  care  of  by  observing  that  the 
result  of  the  product  of  two  terms  with  small  order  variables  is  just  the  linearization  of  the 
final  result  with  respect  to  all  of  the  small  variables.  This  can  be  accomplished  by  taking 
the  Taylor  series  of  the  product  up  to  the  first  order  terms  only.  The  final  result  of  this 
script  is  given  by  equation  (A.  13). 

/*  Construct  the  list  of  error  components  */ 
listl:[ezz,eyz,dxz,exz,dyz,dzz,ezx,eyx,dxx,exx,dyx,dzx]$ 
list2:[ezt,eyt,dxt,ext,dyt,dzt,dxc,dyc,dzc]$ 
list3:[eys,dxs,exs,dys,dzs,dxw,dyw,dzw,ap,ao]$ 
list:append(list  1 ,Hst2,list3)$ 

/*  Declare  Final  Error  Array  */ 
array  (Err,3)$ 

/*  Initialize  all  transformation  arrays  to  0 */ 

array! ioTz.  4,4)$ 

array(izTx,  4,4)$ 

array(ixTt,  4,4)$ 

array! itTc,  4,4)$ 

array(ioTs,  4,4)$ 

array(isTw,  4,4)$ 

array(eoTz,  4,4)$ 

array(ezTx,  4,4)$ 

array(exTt,  4,4)$ 

array(etTc,  4,4)$ 

array(eoTs,  4,4)$ 

array(esTw,  4,4)$ 

for  i:l  thru  4 do 

for  j:  1 thru  4 do 

(ioTz[i  j]  : 0.0, 
izTx[i,j]  : 0.0, 
ixTt[i,j]  : 0.0, 
itTc[i,j]  : 0.0, 
ioTs[i,j]  : 0.0, 
isTw[i,j]  : 0.0, 
eoTz[i,j]  : 0.0, 
ezTx[i,j]  : 0.0, 
exTt[i,j]  : 0.0, 
etTc[i,j]  : 0.0, 
eoTs[i,j]  : 0.0, 
esTw[i,j]  : 0.0)$ 

/*  Initialize  the  ideal  transformation  components  */ 
ioTz[l,l]  : 1$ 
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ioTz[2,2]  : 1$ 
ioTz[3,3]  : 1$ 
ioTz[4,4]  : 1$ 
ioTz[l,4]  : Xl$ 
ioTz[3,4]  : Zl+z$ 
izTx[l,l]  : 1$ 
izTx[2,2]  : 1$ 
izTx[3,3]  : 1$ 

izTx[4,4]  : 1$ 
izTx[l,4]  : X2+x$ 
izTx[3,4]  : Z2$ 
ixTt[  1,1]  : 1$ 
ixTt[2,2]  : 1$ 
ixTt[3,3]  : 1$ 
ixTt[4,4]  : 1$ 
itTc[U]  : 1$ 
itTc[2,2]  : 1$ 
itTc[3,3]  : 1$ 
itTc[4,4]  : 1$ 
itTc[l,4]  : X4$ 
itTc[3,4]  : Z4$ 
i oT  s [ 1 , 1 ] : 1$ 
ioTs[2,2]  : 1$ 
ioTs[3,3]  : 1$ 
ioTs[4,4]  : 1$ 
ioTs[l,4]  : X5$ 
ioTs[3,4]  : Z5$ 
isTw[l,l]  : 1$ 
isTw[2,2]  : 1$ 
isTw[3,3]  : 1$ 
isTw[4,4]  : 1$ 
isTw[l,4]  : xw$ 
isTw[3,4]  : zw$ 

/*  Initialize  Error  Components  */ 

eoTz[U]  : 1$ 

eoTz[l,2]  : -ezz$ 

eoTz[l,3]  : eyz$ 

eoTz[l,4]  : dxz  + ap*dz$ 

eoTz[2,l]  : ezz$ 

eoTz[2,2]  : 1$ 

eoTz[2,3]  : -exz$ 

eoTz[2,4]  : dyz$ 

eoTz[3,l]  : -eyz$ 

eoTz[3,2]  : exz$ 

eoTz[3,3]  : 1$ 

eoTz[3,4]  : dzz$ 
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eoTz[4,4]  : 1$ 
ezTx[l,l]  : 1$ 
ezTx[l,2]  : -ezx$ 
ezTx[l,3]  : eyx$ 
ezTx[l,4]  : dxx$ 
ezTx[2,l]  : ezx$ 
ezTx[2,2] : 1$ 
ezTx[2,3]  •'  -exx$ 
ezTx[2,4]  : dyx$ 
ezTx[3,l]  : -eyx$ 
ezTx[3,2]  : exx$ 
ezTx[3,3]  : 1$ 
ezTx[3,4]  : dzx  + ao*dx$ 
ezTx[4,4]  : 1$ 
exTt[l,l]  : 1$ 
exTt[l,2]  : -ezt$ 
exTt[l,3]  : eyt$ 
exTt[l,4]  : dxt$ 
exTt[2,l]  : ezt$ 
exTt[2,2]  : 1$ 
exTt[2,3]  : -ext$ 
exTt[2,4]  : dyt$ 
exTt[3,l]  : -eyt$ 
exTt[3,2]  : ext$ 
exTt[3,3]  : 1$ 
exTt[3,4]  : dzt$ 
exTt[4,4]  : 1$ 
etTc[  1,1]  : 1$ 
etTc[l,4]  : dxc$ 
etTc[2,2]  : 1$ 
etTc[2,4]  : dyc$ 
etTc[3,3]  : 1$ 
etTc[3,4]  : dzc$ 
etTc[4,4]  : 1$ 
eoTs[l,l]  : 1$ 
eoTs[l,3]  : eys$ 
eoTs[l,4]  : dxs$ 
eoTs[2,2]  : 1$ 
eoTs[2,3]  : -exs$ 
eoTs[2,4]  : dys$ 
eoTs[3,l]  : -eys$ 
eoTs[3,2] : exs$ 
eoTs[3,3]  : 1$ 
eoTs[3,4]  : dzs$ 
eoTs[4,4]  : 1$ 
esTw[l,l] : 1$ 
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esTw[l,4]  : dxw$ 
esTw[2,2]  : 1$ 
esTw[2,4]  : dyw$ 

esTw[3,3]  : 1$ 
esTw[3,4]  : dzw$ 
esTw[4,4]  : 1$ 

/*  Generate  Matrices  */ 
mioTz  : genmatnx(ioTz,4)$ 
mizTx  : genmatrix(izTx,4)$ 
mixTt : genmatrix(ixTt,4)$ 
mitTc  : genmatrix(itTc,4)$ 
mioTs  : genmatrix(ioTs,4)$ 
misTw  : genmatrix(isTw,4)$ 
meoTz  : genmatrix(eoTz,4)$ 
mezTx  : genmatrix(ezTx,4)$ 
mexTt : genmatrix(exTt,4)$ 
metTc  : genmatrix(etTc,4)$ 
meoTs  : genmatrix(eoTs,4)$ 
mesTw  : genmatrix(esTw,4)$ 
oTz  : mioTz. meoTz$ 
zTx  : mizTx. mezTx$ 
xTt : mixTt. mexTt$ 
tTc  : mitTc. metTc$ 
oTs  : mioTs. meoTs$ 
sTw  : misTw. mesTw$ 

/*  Compute  Full  Error  Matrix  */ 

E : (sTwAA-l).(oTsAA-l).oTz.zTx.xTt.tTc$ 

/*  Zero  variables  not  measured  */ 

ezz:0$ 

exz:0$ 

dyz:0$ 

ezx:0$ 

exx:0$ 

dyx:0$ 

ezt:0$ 

eyt:0$ 

dxt:0$ 

ext:0$ 

dyt:0$ 

dzt:0$ 

dxc:0$ 

dyc:0$ 

dzc:0$ 

eys:0$ 

dxs:0$ 

exs:0$ 
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dys:Q$ 

dxw:Q$ 

dyw:0$ 

dzw:0$ 

xw:Xl+X2+X4+x-X5$ 
zw:Zl  +Z2+Z4+z-Z5$ 

/*  Evaluate  E again  */ 
for  i:  1 thru  4 do 

for  j : 1 thru  4 do 

E[i,j]  : ev(E[i,j])$ 

/*  Get  the  linearized  displacement  errors  */ 
for  i:  1 thru  3 do 

Err[i]  : taylor(E[i,4],list,0,l); 
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APPENDIX  B 


This  appendix  lists  the  calibration  coefficients  for  the  original  geometric-thermal  model 
as  well  as  the  adjusted  coefficients.  Since  the  angular  error  equations  are  calibrated  in  arc 
seconds  the  coefficients  are  scaled  to  radians.  Furthermore,  straightness  errors  are 
calibrated  in  micrometers,  and  are  converted  to  millimeters.  All  of  the  other  component 
error  equations  are  calibrated  in  millimeters.  Temperatures  are  in  degrees  Celsius.  This 
discussion  of  units  is  only  given  in  order  to  interpret  the  coefficients  given  below.  The 
methods  described  in  this  report  are  not  dependent  on  these  specific  units. 


B.l  Error  Model  Calibration  Coefficients 

This  appendix  lists  those  initial  coefficients  At]  that  apply  to  most  of  the  machining 

conditions  encountered  on  the  turning  center  portrayed  in  Figure  1.  The  coefficients 
given  here  represent  the  coefficient  values  of  the  component  error  functions  given  in 
Section  3 for  spatial  motion  in  the  negative  x and  positive  z direction.  These  coefficients 
were  developed  during  an  earlier  characterization  of  the  turning  center  (see  [26-27]). 
Spindle  speed  is  assumed  to  be  3000  revolutions  per  minute. 

X Yaw 


An 

A]2 

A|3 

-1.339000E-04 

5.289925E-02 

-3.490320E-01 

Z Yaw 


At] 

A 22 

A23 

1.81 1200E-05 

8.326300E-04 

-6.450000E-01 

X Straightness  with  Z Motion 


A31 

A32 

A33 

A34 

-4.192000E-05 

3.091716E-02 

-4.283536E-01 

1.007243E+01 

41 


X Displacement 


A41 

A42 

A43 

Au 

A45 

3.362300E-07 

2.466 100E-04 

-2.828000E-06 

1.700000E-04 

-1.808700E-03 

Z Displacement 


A51 

A52 

A53 

A 54 

-1.521630E-08 

3.270300E-05 

1.231370E-03 

-3.586880E-02 

Z Straightness  of  X Motion 


A$i 

A63 

^64 

5.097300E-05 

-1.594320E-02 

6.321090E-03 

2.71 1352E-01 

Spindle  Drift 


A71 

A72 

A73 

A74 

-1.020000E-04 

1.20675  IE-02 

-3.894400E-03 

-1.29298QE-Q1 

X Axis  Squareness  with  respect  to  the  Spindle  Axis 


Agi 

A82 

-1.354370E+00 

1.031200E+01 
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Parallelism  of  Z Axis  with  respect  to  the  Spindle  Axis 


A9| 

A92 

A93 

A94 

-7.416000E-01 

3.410200E+01 

4.710000E-02 

-3.93420QE+02 

Machine  Yaw 


Aio.i 

A 10,2 

-5.1 13500E+00 

1.146000E+02 

B.2  Error  Model  Adjusted  Coefficients 


This  appendix  lists  the  updated  coefficients  A . 

X Yaw 


An 

Aj2 

A13 

- 1.5941 62  IE-02 

-6.5335550 

-8.1679140E+02 

Z Yaw 


A21 

A22 

A23 

1.8534552E-03 

5.8945 11 6E-01 

1. 88656 18E+02 

X Straightness  with  Z Motion 


A31 

A32 

A33 

A34 

3.0374870E-04 

-1.4928799E-01 

-3.2036746 

-5.6535790E+01 
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X Displacement 


A41 

A42 

A43 

A44 

A45 

1.3291 142E-06 

4.031 1660E-05 

-1.0762552E-05 

-2.391 838  IE-03 

-6.84 16490E -023 

Z Displacement 


A51 

A52 

A53 

A54 

-5.3957796E-06 

-1.0142461E-03 

-2.1621554E-03 

-1.1814533E-01 

Z Straightness  of  X Motion 


A^i 

A$2 

A^3 

Am 

-3.3490506E-03 

-9.2555965E-01 

-2.7085232 

-8.1 1741 90E+01 

Spindle  Drift 


A71 

A72 

A73 

A74 

-5. 156564  IE-06 

1.487646E-02 

-3.532987  IE-04 

-4.7854670E-02 

X Axis  Squareness  with  respect  to  the  Spindle  Axis 


Agi 

A82 

-1.03561 14E-01 

4.0332 133E+01 
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Parallelism  of  Z Axis  with  respect  to  the  Spindle  Axis 


A91 

A92 

A93 

A94 

-7.6842435E-01 

3.3485040E+01 

-4.9867230E-01 

-4.076 1008E+02 

Machine  Yaw 


Aio.i 

A 10.2 

1.4453392 

2.7857098E+02 
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APPENDIX  C:  Adjusting  for  Too!  Offsets 

All  points  on  a part  are  usually  specified  in  terms  of  part  coordinate  system  points 

(zp,xp).  (C.l) 

Therefore  the  part  origin  must  be  determined  relative  to  the  machine  coordinate  system. 
This  offset  will  be  designated  as: 


(Zoff,Xoff).  (C.2) 


and  is  decomposed  into  two  separate  offsets.  The  first  offset,  (Ztool  set , Xlool  set ) , is  the 

offset  from  the  machine  coordinate  system  to  a point  on  the  turret,  called  the  tool  locating 
point,  to  which  tool  inserts  in  the  turret  are  referenced.  The  second  offset,  (Ztool , Xluol ) , 
is  the  Z and  X values  of  the  tool  point  length  from  the  tool  locating  point.  Thus, 

7—7  + 7 

^ off  ^ tool  _ set  ^ ^ tool' 

lA-  • J ) 

Y — y _l  y 

^ off  tool  _set  ' tool ' 

The  conversion  of  a part  coordinate  to  an  absolute  coordinate  is  then  given  by 


x = xp  +Xoff. 


(C.4) 
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APPENDIX  D:  Determining  the  Machine  Tool  Errors  from  Part 
Measurements 


Although  equation  (18)  associates  the  x and  z residual  machine  errors  at  a general  point 
with  the  perturbation  terms  d , this  equation  cannot  directly  be  used  to  estimate  the 

perturbations.  The  x and  z residual  errors  at  a set  of  sample  points  on  several  parts  must 
be  determined.  These  residual  errors  are  usually  obtained  by  inspecting  parts  on  a CMM 

Assume  that  the  machining  is  done  in,  say,  the  positive  z and  negative  x direction.  This 
is  the  standard  cutting  direction  used  during  this  study.  The  significance  of  this 
assumption  lies  in  the  selection  of  the  characterization  equation  coefficients  given  in 
Appendix  B.l.  Let  there  be  n inspection  points  on  the  surface  of  a part.  One  needs  to 
take  n > 17,  since  there  are  35  unknown  perturbation  terms  and  two  coordinates  per 
inspected  point.  Designate  these  nominal  inspection  points  as 


/ nom  nom 
V 4i 


r nom  nom 


(Dl) 


Suppose  on  a CMM  these  points  are  inspected  and  the  CMM  returns  the  following 
coordinates,  where  the  superscript  p means  inspected, 

(D2> 


Suppose  that  the  process  errors  can  be  identified  at  n inspection  points  on  the  surface  of 
the  part  using  process-intermittent  probing  methods  as  described  in  [17].  If  these  errors 
are  defined  as  dimensional  discrepancies  normal  to  the  surface,  then  the  errors  must  be 
decomposed  into  x and  z components  at  each  of  the  inspection  points.  Thus  we  suppose 
that  we  have  measured  process  error  component  pairs 

(E[p,,EIpi),--,(Kpi,e:pi)(  D3) 


Next,  we  convert  the  inspection  points,  assumed  given  in  part  coordinates  to  absolute 
coordinates  by  computing  using  equation  (C4). 

For  the  absolute  coordinates  (C4)  the  predicted  z and  x errors  are  computed  from 
equation  (2),  using  the  coefficients  from  Appendix  Bl.  Once  this  has  been  done  we  have 
n predicted  machine  error  pairs 


(El,pred,E 


i 

x.pred 


),-,{E%red,Elpred).  (D4) 


The  residual  errors  are  then  computed  at  the  inspection  points  as 
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A E‘ 


A E‘ 


= E‘ 

-E‘ 

~E‘  , 

x,meas 

^X,P1 

x.pred 

- fl 

z,meas 

- E‘ 

^ Z.P 1 

- E‘ 

^ z , pred 

(D5) 


The  subscript  meas  indicates  that  these  errors  are  measured  errors  on  the  CMM  and  the 
subscript  PI  indicates  that  these  errors  are  the  process-intermittent  errors  measured  on  the 
part  while  it  is  on  the  machine  tool. 
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APPENDIX  E:  Numerical  Aspects  of  Pseudo  Inversion 

Equation  (26)  is  in  general  an  overdetermined  system  in  that  there  are  more  rows  than 
columns  or  more  known  parameter  coefficients  than  unknown  parameters.  Ordinary 
inversion  of  7 is  not  possible.  A reasonable  approach  to  obtaining  a “solution”  to  (26)  is 
to  employ  a least  squares  procedure.  If  the  columns  of  7 are  linearly  independent  then 
this  can  usually  be  done  by  fairly  standard  least  squares  software.  However,  if  some  of 
the  columns  are  dependent  then  other  procedures  must  be  applied.  It  is  for  that  reason  we 
introduce  the  idea  of  a pseudo  inverse  that  will  be  described  below. 

In  order  to  “solve”  (26)  for  the  unknown  vector  U we  need  to  examine  some  issues 
involved  with  the  least-squares  solution  of  equation  (26).  A numerical  difficulty  arises  if 
there  are  nonzero  vector  solutions  u of 


Ju=  0 (El) 

There  are  an  infinite  number  of  solutions  of  (El)  since  any  multiple  of  a solution  is  also  a 
solution.  Thus  if  Ui  is  a solution  of  (26)  and  u is  a solution  of  (41)  then,  for  any 
constant  c , the  vector  U x +cu  is  a solution  of  (26)  since  J(Ul  + cu)  = JUl  + cJu  = JUX . 
The  problem  then  is  to  select  from  the  infinite  number  of  solutions  one  particular  “good” 
solution. 

To  analyze  the  difficulty  we  need  to  introduce  some  terminology  from  matrix  theory. 

The  matrix  7 is  a transformation  of  vectors  of  length  35  to  vectors  of  length  2n.  The 
vectors  of  length  35  will  be  called  the  domain  vectors  and  the  vectors  of  length  2n  will  be 
called  the  range  vectors.  The  range  of  J is  a linear  vector  space  in  that  the  sum  of  two 
elements  in  the  range  is  also  in  the  range  as  well  as  a scalar  multiple  of  an  element  in  the 
range  since  7m,  + Juz  = 7(n,  + u2)  and  cJu  - J(cu ) . The  matrix  J can  also  be  written 
in  the  following  form 


j = (j1,j2,-,j3S)  (E2) 


where  7, , for  i = 1,  — ,35 , are  the  columns  of  the  matrix.  Each  of  these  columns  can  also 
be  thought  of  as  a 2n-length  vector.  If  u is  a vector  of  length  35  with  scalar  elements  m( , 
for  / = !,•••, 35  , then  we  can  write,  using  (E2), 


( u ^ 


7m  (7  j,y  2 ’ " " " ’ ^ 35  ) 


1 

M, 


M„ 

35  y 


f>V,  (E3> 


(=1 
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This  shows  that  the  range  of  the  matrix  J is  the  set  of  linear  combinations  of  the  columns 
of  J , thought  of  as  vectors.  A subset  of  the  column  vectors  J.,  for  i = 1,  — ,35 , is  said  to 
be  linearly  independent  if 

0 (E4) 

7=1 


implies  that  u } =0,  for  j = !,•••,&  , where  k < 35 . This  says  that  the  k vectors  (or 

columns)  form  a set  of  coordinate  vectors  for  the  range  space  and  that  k is  the  dimension 
of  the  range  space.  The  number  k is  called  the  column  rank  of  the  matrix  J . If  k is  35 
then  the  matrix  is  said  to  be  of  full  rank.  The  null  space  of  the  matrix  J is  the  set  of  all 
vectors  u that  satisfy  Ju  = 0.  This  set  is  also  a vector  space.  Its  dimension,  also  called 
its  nullity , is  the  largest  set  of  linearly  independent  vectors  that  satisfy  Ju  = 0.  A 
fundamental  result  in  matrix  theory  is  that  the  sum  of  the  rank  and  the  nullity  of  a matrix 
is  the  dimension  of  the  domain  space. 

In  order  to  determine  a least-squares  solution  of  equation  (26)  some  measure  of  the 
difference  between  the  right  and  left-hand  sides  of  (26)  must  be  defined.  The  usual 
measure  or  metric  used  extends  the  idea  of  Euclidean  distance  to,  in  this  case,  vectors  of 
length  2 n,  where  n is  the  number  of  probed  points  on  the  workpiece.  We  will  denote  the 
length  or  norm  of  a 2«-vector  by 


H = <ES) 

The  least-squares  problem  associate  with  equation  (26)  is  to  find  the  2/i-vector  that 
minimizes  the  norm  of  the  difference  between  the  right  and  left-hand  sides  of  (26).  This 
is  symbolically  written  as 

mini  Ju  - v||  (E6) 

u 1 


where  v is  a right-hand  side  vector  for  (26).  From  linear  regression  analysis  the  solution 
of  (E6)  depends  on  solving  the  system 

JTJu  = JTv  (E7) 

where  the  superscript  T represents  the  transpose  of  the  matrix.  There  are  two  cases  to 
consider:  (1)  rank(J)  = 35  and  (2)  rank(J)  = k < 35.  If  the  rank  of  J is  35,  then  the  rank 

of  JT  J is  35  and  it  is  a square  matrix.  Therefore  JT  J has  an  inverse  and  (El)  is 
uniquely  solvable.  This  is  the  traditional  least-squares  solution  of  (E6).  However,  if 
rank(J)  = k < 35  , then  the  rank  of  JT  J is  k(<  35),  which  means  that  the  nullity  of  JT  J 
is  nonzero  and  (E7)  has  an  infinite  number  of  solutions  or  an  infinite  number  of  vectors 
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that  satisfy  (E6).  The  object  then  is  to  select  an  appropriate  solution  for  (38)  and 
therefore  ultimately  for  (26). 

A significant  method  for  decomposing  matrices,  such  as  J,  is  called  the  Singular  Value 
Decomposition  and  is  strongly  related  to  the  eigenvalue-eigenvector  decomposition  of 
J 7 J . The  singular  value  decomposition  of  the  2 n by  35  matrix  J of  rank  k is  given  by 

J = ODWt  (E8) 

where  O,  a 2 n by  In  matrix  and  W,  a 35  by  35  matrix,  are  orthogonal , which  means  they 
satisfy 


00T  -Ot O = / 
WWT = WTW  =/ 


(E9) 


and  D is  a 2 n by  35  diagonal  matrix  with  nonnegative  diagonal  elements  arranged  to  be 
nonincreasing  from  upper  left  to  lower  right.  Some  of  them  in  the  lower  right  can  be  0. 
The  relation  to  the  matrix  J7  J is  given  by 

JtJ=WDtDWt  (E10) 

which  is  not  hard  to  show  by  multiplying  the  transpose  of  (E8)  times  (E8)  and  using  (E9). 
This  shows  that  the  singular  values  of  J are  the  positive  square  roots  of  the  eigenvalues  of 
JTJ  . 


The  singular  value  decomposition  of  J can  be  used  to  solve  the  least-squares 
minimization  problem  (E6)  by  using  the  following  result  (See  [24]). 

Suppose  that  J is  decomposed,  using  the  singular  value  decomposition  to  the  form  (E10) 
where  the  matrix  D can  be  written  in  the  block  form 


D = 


(D\\ 

0 


0" 


(Ell) 


where  the  submatrix  Dn  is  a diagonal  kxk  matrix  with  the  nonnegative  diagonal 
elements,  containing  the  singular  values,  arranged  to  be  nonincreasing  from  upper  left  to 
lower  right.  Introduce  two  new  variables  g and  y by 


0Tv  = g = 
WTu  = y = 


V#2y 


(E12) 
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where  g, , y,  are  vectors  of  length  k,  the  rank  of  Dn , g2  is  a vector  of  length  2 n-k  and 
y2  is  a vector  of  length  35 -L  Note  that  at  this  point,  since  u is  unknown,  y is  unknown. 
Now  let  5?i  be  the  unique  solution  of  the  matrix  equation 

Dnyi=8i  CE13) 

Then  the  first  conclusion  is  that  all  solutions  of  the  minimization  problem  (E6)  are  of  the 
form 


u=W 


(El  4) 


where  y2  is  an  arbitrary  vector  of  length  35-£.  Next,  any  vector  of  the  form  (E14) 
generates  the  same  residual  vector  r given  by 


r = v - Ju  = O 


0\ 

\8lJ 


(El  5) 


This  says  also  that  the  vector  of  the  difference  between  any  two  vectors  of  the  form  (E14) 
lies  in  the  null  space  of  J.  Therefore  any  two  vectors  of  the  form  (E14)  differ  by  a vector 
in  the  null  space  of  J.  The  norm  of  r satisfies  ||r||  = ||v  - Ju\  = ||g2||  and,  finally,  the  unique 
solution  of  the  minimization  problem  (38)  with  the  minimum  length  is  given  by 


u = W 


V J 


(El  6) 


The  fact  that  (E16)  is  the  unique  solution  of  minimum  length  of  (E6)  permits  the 
definition  of  a concept  of  generalized  inverse  of  the  matrix  J even  though  it  is  not  a 
square  matrix.  In  particular  the  following  result  also  holds.  If  the  matrix  J is  written  in 
the  form  of  its  singular  value  decomposition  (E8)  where  the  diagonal  matrix  D is  written 
in  the  block  form  (Ell)  then  the  unique  minimum  length  solution  of  the  minimization 
problem  (E6)  can  also  be  written  as 


u =W 


fD~ 1 0^ 


V 


11  v Otv  (E17) 

0 oj 


by  combining  (El 2),  (E13)  and  (E16).  Now  define 


( r\- 1 


i+  = w 


o 

0 0 


'v 


0T  (E18) 
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The  matrix  J + is  called  the  generalized  inverse  of  J. 

Although  the  matrix  Du  has  all  of  its  diagonal  elements  nonzero  so  that  it  has  the  form 


A,  = 


(dx 

0 

0 


0 "l 

0 

dk 


(E19) 


care  must  be  taken  when  the  inverse  is  computed  numerically.  Theoretically  the  inverse 
of  Dn  has  elements  of  the  form  — ,i  = !,•••,&,  but  the  values  of  d may  be  sufficiently 

d, 

small  so  that  their  reciprocals  are  numerically  meaningless.  Thus,  when  the  generalized 
inverse  of  a matrix  is  computed  the  diagonal  matrix  of  singular  values  Du  is  usually 
redefined  to  be  the  matrix  that  contains  only  the  singular  values  whose  absolute  values 
are  larger  than  some  tolerance.  All  others  are  set  to  zero,  so  that  from  a numerical  point 
of  view  the  rank  of  Du  may  be  less  than  k.  A tolerance  that  is  used  often  in  determining 
those  singular  values  that  are  to  be  eliminated  is  given  by 

tol  = 2n||y||f  (E20) 


where  the  norm  of  a matrix  is  computed  in  a similar  manner  to  the  norm  of  a vector  in 
that  the  square  root  of  the  sum  of  the  squares  of  all  of  the  elements  is  taken  as  the  norm 
and  e is  called  the  machine  epsilon  and  is  that  positive  finite  machine  precision  number 
that  satisfies  1 + e = 1 . For  an  ordinary  PC  using  double  precision  the  number  is 
approximately  2.2E-16. 
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